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ABSTRACT 

We simulate one-dimensional arrays of tunnel junctions using the kinetic Monte Carlo method to study charge filling behaviour 
in the large charging energy limit. By applying a small fixed voltage bias and varying the offset voltage, we investigate this 
behaviour in clean and disordered arrays (both weak and strong disorder effects). The offset voltage dependent modulation of 
the current is highly sensitive to background charge disorder and exhibits substantial variation depending on the strength of the 
disorder. We show that while small fractional charge filling factors are likely to be washed out in experimental devices due to 
strong background charge disorder, larger factors may be observable. 


Introduction 

Since the prediction 1-3 and subsequent observation 4,5 of single charge tunnelling in small metallic tunnel junctions, the field of 
low temperature nanoelectronics has developed rapidly. Fabrication of ultrasmall metallic islands allows the discrete nature 
of the island charge states to be observed. This leads to the suppression of conduction at low bias voltages, due to Coulomb 
blockade. However, introducing multiple junctions increases the complexity of the charge dynamics significantly. Exploring 
and understanding the structure and dynamics of low energy excitations in multi-junction systems and the effect of random 
offset charges has become important in improving the accuracy and design of these circuits. 

An example of such a multi-junction circuit is a tunnel junction array which consists of a chain of small metallic islands 
separated by tunnel barriers, see Fig. 1. Transport through the array can be correlated due to the interplay between the bias 
voltage across the array and the Coulomb repulsive interaction between charges. This interaction decays exponentially over a 
characteristic length A that is determined by the ratio of the circuit capacitances. Consequently, in the low bias regime, charges 
within a junction array tend to be equidistant and in long arrays, at the onset of charge injection, a charge pattern quickly 
emerges. 

These patterns are dependent on the charge filling within the array, which is an additional modifier of the current. Different 
patterns emerge when certain filling factors are enforced by applying an offset voltage U , which changes the chemical potential, 
resulting in variation in the charge periodicity within the array. For example, a specific charge configuration within the array 
can mean the difference between a maximal current signal and current quenching. 

Owing to their nanoscale geometry and as a result, large charging energy, these systems are highly sensitive to even small 
fractional random offset charges in the surrounding substrate. The effects of random charge disorder have been investigated in 
various single charge devices including transistors 6, and arrays, 8- however, these studies focussed primarily on the shift in 
the threshold voltage and/or properties of the device noise. 

The effect of background charge polarisation in the substrate of a tunnel junction array was investigated in Ref. 12 and the 
effect of this background was found to be non-trivial. That work involved putting a random fractional charge on a random 
number of islands in the array and allowing it to come to equilibrium at zero bias. However, our work concerns the transport 
properties at small fixed bias. Experiments thus far in this regime are consistent with maximal disorder behaviour. 13 In 
sufficiently long devices maximal disorder is expected to completely suppress the U -dependence of the transport properties by 
self-averaging. The devices fabricated so far however are too short for full self-averaging as seen in standard condensed matter 
systems. It is therefore important to study the ^/-dependence in short arrays even in the case of maximal disorder. 

Improvements in fabrication technologies may reduce disorder in these devices. Therefore, it is important to understand the 
effects of varying levels of disorder and to which degree disorder would have to be reduced to see real improvements in the 
behaviour of these devices. We investigate charge filling factors by simulating a junction array model in the charging energy 
limit (negligible Josephson energy) with zero disorder. We then extend this model to include various strengths of background 
charge disorder and consider the effect disorder has on the charge filling structure. We specifically focus on the signatures of 



Figure 1 . Top: Schematic of the junction circuit we simulate, consisting of a linear chain of N +1 tunnel junctions each with 
junction capacitance Cj. A symmetric bias V is applied across the array with an additional offset U relative to ground and each 
island is capacitively coupled to a ground plane through C G - Bottom: I-V-U characteristics for 20 different offset voltages from 
U = 0 to U = 0.5 e/Cc- Due to the voltage offset U , the I-V responses are not simply linear within this V range and the 
threshold voltage is also [/-dependent. The current at fixed offset voltage bias is also plotted (magenta line), for V = 0.1 e/Cc . 


filling factors and background disorder in the current-voltage characteristics of the device. 


Clean array 

We begin by considering an ideal model wherein the array is clean and homogeneous - there is no background charge disorder. 
We consider a one-dimensional linear array consisting of N = 50 islands each with a junction capacitance Cj = 50 aF, an 
effective capacitance to ground C G — 2 aF and an effective electron temperature of T e = 30 mK. The array is driven by a 
symmetric voltage bias consisting of a small fixed bias V and varying offset voltage U, see Fig. 2. The current-voltage 
characteristics are non-monotonic functions of U due to charge filling factors within the array. We consider the normal 
conducting regime where the current is comprised entirely of electrons due to single electron tunnelling. The charge interaction 
length can be approximated by A = yjC j/Cq , where A = 5 (or A = 15, where specified) for our simulations. The Hamiltonian 
consists of a charge interaction term and two source terms and is given by 
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for a particular charge configuration q. The two source terms 5i and <5/v are equal to one for site 1 and A, respectively and zero 
everywhere else. The electrostatic interactions of the array are described by the N x N capacitance matrix 
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Figure 2. Current response as a function of U and V in steps of 0.025 e/Cc. The fine structure is sharpest at approximately 
V = 0.1 e/Co- While the finer charge filling structure is lost at large V, this structure is preserved for smaller V before 
becoming washed out by noise. The period of the current in U is given by e/Cc (dashed-dotted line). 


and the charges can incoherently tunnel from island to island parameterised by a tunnelling resistance R t . We simulate the 
charge dynamics of the array using the kinetic Monte Carlo (KMC) method. 3,14-17 At each Monte Carlo step, all possible 
single electron movements are computed. The rate for each of these transitions is given by 18 

r(5£) = Jr, Cl dEf{ - E)[l -f( E+5E )} 

where 8E is the energy difference between initial and final charge configurations and f(E) is the Fermi function. In general, 
the influence of the circuit impedance on the transition rate can be included via P(E)-theory. 18 However, as Cg/Cj 1, we can 
consider the junction array to be a low impedance environment. 19 In this limit, P(E) ~ 8(E), resulting in the expression given 
in Eq. 3 for the transition rates. These rates are used to select which transition occurs, based on a weighted probability. The 
total current is then computed based on the net movement, left and right, of charges over many Monte Carlo steps. To ensure 
the array has equilibrated, Monte Carlo steps of 10 5 -10 6 are used before collecting statistics with an additional 10 5 -10 6 steps 
and ensemble averaging over multiple data sets when considering disorder. For a detailed review of the KMC algorithm, we 
refer the reader to Ref. 20 or Ref. 21 . 

We can also calculate the analytical threshold voltage for a current to flow through the array (at T e = 0). 16,22 In the limit of 
long interaction length (A 1) and assuming an initially empty array 

= 2\/CjCc: (4) 

where e is the elementary charge. The threshold voltage and the current response are highly variable, depending on the 
characteristic value of the offset voltage U which sets the filling fractions within the array, as shown in Fig. 1. This behaviour is 
periodic in U, which is the focus of this paper, and is shown for half a period via the magenta line in Fig. 1. 

The {/-dependence can also be seen in Fig. 2 where we plot the current response as a function of U and V. The current is a 
periodic function of the offset voltage U for all simulated values of V, where the period of the current in U oc e/Cc (the effective 
voltage on the ground capacitor). Therefore, the offset voltage required for an integer filled array (point of zero current) is 
filled = Qe/CG, where Q G Z. Due to particle-hole symmetry, the response is symmetric about the point U = e/2Cc • The 
fine structure, seen at smaller voltages, is lost as V is increased and the response flattens out and becomes more rounded. It is 
this fine structure that is a signature of charge filling fractions. 

In order to identify filling factors (p/q), we compute the time averaged charge-charge autocorrelation function (fi/(0)fii+ m (0)) 
to measure the statistical correlations within the charge distribution. As the length of the array is not infinite, boundary effects 
are significant and the periodicity of the filling factors is strongest at the centre of the array. 

From the current modulation in U, we can identify charge correlation patterns, i.e., peaks, troughs and dislocations in the 
patterns, as shown in Fig. 3. We can understand this structure by considering the filling factors of the array. Current blockade 
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Figure 3. I-U response at A = 5 and A = 15 for V = 0.1 e/Cc and V = 0.0226 e/Cc , respectively. Dominant p/q filling 
fractions are labelled for specific features. At A = 5, mixing of different filling fractions is seen due to the short interaction 
length-scale. The charge-charge autocorrelation functions (tii(0)ni+ m (0)) , shown for the first 15 sites, provide a quantitative 
approach to identifying the filling fractions within the array. The autocorrelation shows the washing out of the filling factors as 
they approach 1/A = 1/5. The colour bar denotes the magnitude of the autocorrelation function. 


occurs at integer filling, i.e., the array is uniformly filled and it is difficult to move this pattern forward by injecting another 
charge or otherwise disrupting the pattern. In other words, integer filling factors correspond to states in which the charge 
configuration is rigid, therefore, the system requires significantly more energy to break this pattern and allow conduction to 
commence. 

Fractional filling, on the other hand, allows current to flow more easily due to the presence of unoccupied sites between 
current carrying charges. For a perfectly periodic charge distribution that matches the boundary conditions of the array, we also 
see quenching. However, when the periodicity of the charges is disrupted, either via mismatch with the boundary or (near) 
degeneracy of two different charge configurations, current can flow more easily via interconversion of these charge states. As 
an example, periodic p/q filling fractions (e.g., 1/3 filling) are more difficult to move forward, resulting in current quenching. 
The combination of two different filling fractions (e.g., 1/3 and 1/2 filling) is easier to transfer through the array because the 
two factors can be interchanged. Therefore, as a function of U, we see current maxima and minima as aperiodic and periodic 
charge patterns, respectively. In Fig. 3, this effect can clearly be seen for the case of A = 15. 

However, when the average spacing between charges approaches the interaction length A, we see a greater mixing of the 
filling fractions and the troughs in the current magnitude shift. As p/q—> 1/A, we see peaks rather than troughs at rational 
values of p/q. This can be more clearly seen in Fig. 3 for the case of A = 5. When the mean charge separation is < 1/A, for 
example, 1/2, mixing of the charge configurations is weaker. Mixing at small A = 5 can also be seen in the autocorrelation 
functions in Fig. 3, where the fractions are no longer clear and become blurred, in contrast with the clarity of the same filling 
factors in the autocorrelation for large A = 15. 

We observe splittings associated with defects in the dominant patterns, for instance, at A = 5, 1/2 filling in Fig. 3. The way 
in which these patterns match against the boundary conditions splits the equivalency of the different degenerate configurations. 
For example, at A = 5, 1/2 filling, we see that slightly above or below p/q — 1/2 gives a slightly lower current (the trough- 
peak-trough structure). This feature is not seen atA = 15, p/q = 1/2 because at large interaction lengths averaging over a 
greater number of sites washes out this effect. 

This is controlled by how easy it is to break the pattern, i.e., by adding an additional electron (or lack thereof) to make the 
pattern shift by half a period. As an example, the (metastable) state q T = 10101011010101 is slightly harder to move along the 
array than the ‘perfect’ pattern of q T = 10101010101010. Therefore, depending on the matching of the boundary conditions 
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Figure 4. I-U characteristics at V = 0.1 e/Cc as a function of site disorder strength in steps of 0.01 from rj = 0 to rj = 0.1. 
Each plot is ensemble averaged over 50 different disorder realisations. Even at very low disorder (rj = 0.01), while the 
structure remains largely unchanged, the current magnitude is reduced. The magnitude of the current continues to decrease 
with increasing disorder strength. 


and the presence of defects changes whether we see a peak or a trough. 

Effects of disorder 

Experimentally, imperfections or a degree of disorder is always prevalent in the system (e.g., in the substrate and within the 
junctions themselves), therefore, the inclusion of disorder makes our analysis more applicable to experimental observations. 
Furthermore, an important question is which features of this analysis are preserved in an inhomogeneous array. In this section, 
we extend our model by including various background charge disorder and investigating the weak and maximal disorder limits. 
We simulate boxed disorder (a uniform distribution) by adding a random fractional charge to each site, q G [— r\e, +r}e]. The 
maximum/minimum value of disorder on any given site is given by ±r/e, where T] indicates the width of the disorder distribution. 
This type of potential disorder simulates static background charges and is often used to study mesoscopic systems. 9,12,23 

The effect of the background on the current-voltage response can be seen even at nominal disorder (e.g., q = 0.01). While 
the structure remains largely unchanged, the current magnitude steadily decreases, as shown in Fig. 4. As the magnitude of 
static offset charges increases, the response becomes noisier and eventually only the most prominent structure remains (e.g., 
integer and 1/2 integer filling). However, even at very small current magnitude (strong disorder) periodic structure in U remains. 
The current scales as rj because disorder constricts the conduction channel, making it more difficult for charge to flow (as 
opposed to a clean array). At strong disorder, a larger V is required to overcome the random static background charges and 
force current through the array. 

Maximal disorder should be reached at T] = 0.5 as this value of disorder corresponds to the point where any value larger 
than 0.5 can be negated by one additional charge tunnelling onto the island. Indeed, at 7] « 0.5, our model shows maximal 
disorder and all of the structure is lost and charge filling factors are no longer seen. In Fig. 5, we show that the current is 
approximately constant across U at rj =0.5, indicating that maximal disorder has been reached (for this value of V ) and a 
lack of {/-dependence. Then at rj = 0.6 — 0.9, the current oscillates again, with opposite curvature (on average) to the T] < 0.4 
currents, with the current once again becoming approximately constant at rj = 1. In order to understand this behaviour, we 
consider three different values of disorder; less than, equal to and greater than maximal disorder, depicted in Fig. 6. When the 
disorder distribution is broader than maximal disorder (rj = 0.5), adding an individual whole charge causes the distribution 
to wrap back around (or fold back) into the interval ±0.5. This folded disorder distribution can be mapped back to the 
original box-distribution form with an additional ‘base’ contribution spanning the whole [—0.5,0.5] interval by a shift of 0.5e 
(corresponding to the opposite current curvature observed in Fig. 5). 

In Fig. 7, we plot the current response as a function of T] for three different values of A at both the maximal and minimal 
I-U points. The applied voltage is scaled such that V = 0.8V t h for each value of A, where Vth is calculated via Eq. 4. In the 
clean array limit, the maximum current depends on A because the current response in the transport regime is not exactly Ohmic, 
but shows a linear current-voltage dependence with a non-zero current offset. The maximum/minimum currents as a function of 
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Figure 5. I-U response at V = 0.2 e/Cc for various disorder strengths in steps of 0.1 from rj = 0 to rj = 0.6 ,17 = 1 is also 
shown (dotted line). Each plot is ensemble averaged over 50 different disorder realisations. Points of maximal and minimal 
current at zero disorder are shown (dashed lines). The system reaches maximal disorder at 77 =0.5. When 0.5 < rj < 1, the 
current exhibits oscillatory behaviour, however, the system reaches maximal disorder again at 77 = 1. The mean current across 
U and the variance in the disorder realisations for 0.1 < 77 < 1 is shown via error bars which depict one standard deviation. 



Figure 6. Pictorial representation of the probability distribution of the disorder as a function of 17 showing maximal disorder 
and the folding of the charge disorder back in on itself to maximal disorder (77 = 0.5). For very broad distributions, 17 >> 0.5, 
the 77 -dependence will abate due to decreasing ‘base’ to ‘top’ ratio. 


U for all studied values of A equilibrate at approximately rj = 0.5, i.e., maximal disorder. This is because disorder suppresses 
the current variations as a function of U, seen in Fig. 5. The oscillatory behaviour of the current at 0.5 < 77 < 1 can also be 
seen, with the currents at 77 =0.5 and 77 = 1 equal for each value of A due to maximal disorder folding. Note that throughout 
this discussion A < N. 

To better understand the role of the background charge disorder, we introduce the concept of a disorder correlation length, 
the distance over which charges on adjacent sites are still correlated. Beyond this length, the background disorder has the effect 
of randomising the electrostatic environment that the charges see and therefore, islands separated by more than this length can 
be considered independent. Assuming that the cumulative effect of charge disorder can be thought of as a random walk style 
process , 24 we define the disorder correlation length as 


V V^con — 2 


(5) 


where 77 is the value of charge disorder. When rj >0.5, L corr < 1, which corresponds to a delta-correlated disorder since the 
minimal length-scale of the system is the inter-site-distance L = 1. Such a delta-correlated disorder corresponds to the maximal 
disorder model . 13,24,25 Conversely, 77 —^ 0 corresponds to L corr —» oo, i.e., the clean array limit. 


6/9 











































^corr 


oo 25 6.25 2.778 1.563 1 0.694 0.510 0.391 0.309 0.25 



Figure 7. Current response as a function of 77 for A = 5,10,15, ensemble averaged over 100 different disorder realisations. 
The applied voltage is scaled such that V = 0.8V t h for each value of A. The current response is shown for both 7 max (triangles) 
and /min (circles), i.e., the points of maximal and minimal (zero) current, respectively, in the I-U response. The values of L corr 
corresponding to the three values of A (dashed lines) and maximal disorder (solid line) are also shown. 


When V is constant, L corr = 20,35,50 have near identical I-U characteristics as a function of N (data not shown). In contrast 
to the A-dependence of the threshold voltage, 26 the system exhibits similar behaviour whether L corr is equal to, less or greater 
than A. At A = L co lix , the onset of transport changes from a boundary to a bulk effect and consequently, the dependence of V t h 
on A changes. Here, however, we are considering a tunnelling array that is already in the transport regime and the modulation 
of the current due to a change in filling factors is naturally a bulk effect. 

Discussion 

In conclusion, we have studied charge filling factors as they manifest in a (unrealistic) clean tunnel junction array. These effects 
are predicted to be difficult to observe experimentally as the intrinsic disorder in experimental devices likely approaches or 
exceeds the limit for their observation. However, while small fractional filling factors are expected to be unobservable in present 
day devices, our model predicts that observation and reproducibility of larger factors may be possible (excluding maximally 
disordered arrays). For example, current blockade near integer filling should be quite robust to static offset charges. Signatures 
of charge filling factors could also be observable in other one-dimensional systems. The clean array analysis is applicable to 
arrays of quantum phase slip elements 27 where the effects of localised flux disorder should be much weaker and should display 
longer disorder correlation lengths. The maximal disorder analysis is more applicable to Josephson junction arrays, where 
background charge offsets dominate the behaviour. 
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